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ABSTRACT 

Finding electromagnetic (EM) counterparts of future gravitational wave (GW) sources 
would bring rich scientific benefits. A promising possibility, in the case of the coales- 
cence of a super-massive black hole binary (SMBHB), is that prompt emission from 
merger-induced disturbances in a supersonic circumbinary disk may be detectable. We 
follow the post-merger evolution of a thin, zero- viscosity circumbinary gas disk with 
two-dimensional simulations, using the hydrodynamic code FLASH. We analyze per- 
turbations arising from the 530 km s _1 recoil of a 10 6 M Q binary, oriented in the plane 
of the disk, assuming either an adiabatic or a pseudo-isothermal equation of state for 
the gas. We find that a single-armed spiral shock wave forms and propagates outward, 
sweeping up ~ 20% of the mass of the disk. The morphology and evolution of the 
perturbations agrees well with those of caustics predicted to occur in a collisionlcss 
disk. Assuming that the disk radiates nearly instantaneously to maintain a constant 
temperature, we estimate the amount of dissipation and corresponding post-merger 
light curve. The luminosity rises steadily on the time-scale of months, and reaches 
few xlO 43 erg/s, corresponding to « 10% of the Eddington luminosity of the central 
SMBHB. We also analyze the case in which gravitational wave emission results in a 
5% mass loss in the merger remnant. The mass-loss reduces the shock ovcrdcnsitics 
and the overall luminosity of the disk by w 15 — 20%, without any other major effects 
on the spiral shock pattern. 
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1 INTRODUCTION 

The gravitational waves (GWs) produced during the late 
stages of the merger between super-massive black holes 
(SMBHs), with masses of ~ (10 4 -10 7 ) M /(l + z), out to 
redshifts beyond z « 10, are expected to be detectable in 
the next decade by the Laser Interferometric Space Antenna 
(LISA) satellite. While the GW signatures themselves will 
be a rich source of information, identifying the electromag- 
netic (EM) counterpart of the LISA source would open up a 
whole new range of scientific opportunities, from black hole 
astrophysics to fundamental aspects of gravitational physics 
and cosmolog y (|Kocsis et al.|[2007l |Haiman et "aLl|2009a| 
Phinney|2009) |Bloom et al.||2009[). 



Whether the EM counterparts of the SMBHBs detected 
by LISA can be uniquely identified depends on the accuracy 
of localization by LISA, and on the nature of the EM emis- 
sion. For the typical SMBHB detected by LISA at z « 1 - 2, 
the sky position and the redshift will be known to the accu- 
racy of S(AQ) ~ few x 0.1 square degrees, and 5z « few x 
0.01, respectively (the latter limited by weak lensing errors; 
e.g. |Holz fc Hughes||2005| |Kocsis et aT1|2006[). Within this 



three-dimensional error volume, there will be of order ~ 100 
candidate galaxies, at the optical magnitude limit expected 
to host such SMBHBs ( Kocsis et al.|2008 l. If the coalescing 
SMBHs themselves produce bright emission comparable to 
luminous quasars, or are associated with some other, sim- 
ilarly rare subset comprising < 1% of all galaxies (such as 
ultra-luminous infrared galaxies), then a unique counter- 



part may be identified among these candidates (Kocsis et 



al. 20061. However, having a prediction for the spectrum 
and the light-curve of a coalescing SMBHB - and therefore 
knowing what characteristic signatures to look for - would 
make such identifications both more likely and more reliable. 

A promising possibility is that the coalescing SMBHB 
produces a variable or transient EM signal, which can be 
uncovered by suitably designed EM observations, either con- 



currently with or following the LISA detection ( Kocsis et al. 
|2008}|Lang fc Hughes|2008l|Haiman et al.|2009a[ ). The dense 
nuclear gas around the BH binary is expected to cool rapidly, 
and settle into a rotationally supported, thin circumbinary 
disk (e.g. |Barnes||2002] |Escala et aT1[2u05| . The evolution 
of a SMBHB embedded in such a disk has been studied in 
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various idealized configurations (e.g. 



Armitage fc Natarajan| 



2002| Liu, Wu fc Cao 2003; Milosav 



Dotti et al. 2006; MacFad yen 
20091 ICuadra et al.|2009||Haiman et al.|2009b| 



jivic fc Phinney||2005[ 
Mil osavljevic|2008||Hayasaki| 
Generically, 



at small orbital separations when the binary is detectable 
by LISA, the orbit decays rapidly due to GW emission. If 
the disk is thin, the binary torques create a central cavity, 
nearly devoid of gas, within a region about twice the or- 
bital separation (e.g. Artymowicz fc Lubowj 19*94 1. Whether 
the decaying SMBHBs produce bright emission during this 
stage is not well understood. If the central cavity were truly 
empty, no gas would reach the SMBHBs, and any emission 
produced farther out in the disk would likely be weak. On 
the other hand, numerical simulations suggest residual gas 



inflow into the cavity ( Artymowicz fc Lubow| 1996| |Mac- 



Fadyen fc Milosavljevic||2008| |Hayasaki et aL|[2007| [2008 



Cuadra et al. 20091, which may plausibly accrete onto the 



BHs, producing non-negligible EM emission. 

A different possibility, and the focus of the present pa- 
per, is that variable EM signatures are produced in the 
gas disk promptly after the coalescence of the SMBHB. The 
burst of GWs emitted during the last stages of the coales- 
cence results in a corresponding, nearly instantaneous re- 
duction in the binary's rest mass. Furthermore, when com- 
pact objects coalesce asymmetrically, the linear momentum 
carried away by the GWs imparts a "kick" to the center of 
mass of the system ( Bekenstein|1973 Fitchett|1983 1. Recent 
break-through in numerical relativity has allowed accurate 
calculations of these effects, showing that the mass loss is 
typically several percent (e.g. Tichy fc Marronetti|2008" and 
references therein), while kick velocities are typically sev- 
eral hundred km s _1 , but can be as high as 4,000 km s _1 
( Baker et al.||2006| |2007| [2008] [Campanelli et~aT||2007a|b[ 



Gonzalez et al.||2007a|b| |Herrman et al.||2007a|b|c| |Koppitz| 



et al. 2007 1. The circumbinary gas will respond promptly (on 



the local orbital timescale) to such dynamical disturbances. 
Importantly, the gaseous disk outside the inner cavity is ex- 
pected to cool efficiently and become geometrically thin, im- 
plying that the orbital motion of the gas is supersonic. This 
gas is therefore susceptible to prompt shocks, which could, 
in principle, produce a detectable transient EM signature 
( jMilosavljevic fc Phinney|2005| . 



Lippai et al. (2008) considered the motion of collision- 



less test particles in such a disturbed disk. As long as the 
particles remain bound to the central SMBHB, they fol- 
low elliptical Kepler orbits (in the inertial frame of the 
SMBHB). These orbits cross, and produce a characteris- 
tic, outward-propagating spiral caustic. While these results 
were obtained in a pressureless "dark matter" disk, they 
suggest that shocks, with a similar pattern, will indeed arise 
in the gas, on times-scales of ~weeks to ~months after the 
coalescence. Similar conclusions were reached by Shields & 
Bonningl (120081) and by ISchnittman & Krolikl (|2008|), using 



N-body simulations and semi-analytical arguments, respec- 
tively, to show that a bright, post-merger "flare" may oc- 
cur. Both of these studies focused on the evolution of disks 
around more massive BHs on longer (~ 10 4 yr) time-scales, 
and proposed detecting the flare by monitoring a population 
of active galactic nuclei (AGN). By comparison, our anal- 
ysis here focuses on the disks around lower-mass BHs and 
on the shorter time-scales of several weeks to a year, which 



are relevant for follow-up observations triggered by LISA 
detections. 

Motivated by the above, in the present paper, we fol- 
low up on earlier work, and compute the response of the 
gas disk including the effects of gas pressure using hydrody- 
namical simulations. Our main goals are (i) to assess shock 
formation in hydrodynamical disks, and (ii) to estimate the 
amount of dissipation and the resulting light-curve pro- 
duced by the disk. In particular, for the latter, we will use 
simulations with a pseudo-isothermal equation of state. This 
corresponds implicitly to strong dissipation, and should rep- 
resent an approximate upper limit on the disk luminosity. 
Our computations are performed with the publicly available 
code FLASH ( |Fryxell et al.|2000| ). 

Two other, independent studies have recently used 3- 
dimensional simulations to address the response of gas disks 



to mass-loss (O'Neill et al. 20091 and to both mass-loss 
and kicks ( |Megevand et a l. 2009). The main difference be- 
tween our analysis and these previous works is our inclusion 
of runs with a pseudo-isothermal equation of state, which 
modifies, qualitatively, the expected EM signature. In par- 



ticular, O'Neill et al. (2009) and Megevand et al. ( 2009 ) both 



run simulations with an adiabatic equation of state, and find 
that in most cases, the gas disk typically dims following the 
merger. We observe a similar effect in our adiabatic runs, 
and attribute it to an overall dilution of the gas density. 
However, the gas develops significantly larger density con- 
trasts in our pseudo-isothermal runs, and we argue that this 
can result in a significant brightening of the system. Another 
important difference between our study and those of |Q'Neil"l| 
et al. (2009) and Megevand et al. (20091 is that we perform 



two-dimensional simulations, whereas O'Neill et al. (20091 



and |Megevand et al?] ( |2009| | both utilized 3D simulations 
While we have sacrificed resolving the vertical disk struc- 
ture, this allows us to simulate a disk that extends to much 
larger radii (10 4 Schwarzschild radii), and follow the entire 
disk for a much longer period (~ 1 yr for a 10 6 Mq bi- 
nary), in order to cover the regime relevant for follow-up 
observations of LISA events. Given the expected size of the 
inner cavity in the disk, ~ 100-R S , studying these outer re- 
gions of the disk, and following the disk evolution on the 
correspondingly longer time-scale, is especially important. 
Further differences between our study and those of |Q'NeiIT| 
et al.| ( 2009 1 and Megevand et al. ( 2009 ) will be discussed in 



3761 below. 

The rest of this paper is organized as follows. In § [2] 
we describe our computational setup, including details of 
the simulations, and initial conditions. In § [3] we present 
and discuss our main results. These include the impact of 



the kick (§ 3.1| and the mass-loss (as well as both effects 
in combination; § 3.2 l, using constant surface density disks. 



In § |3.3| we discuss our estimate of the post-merger light 
curves. In § |3.4| we repeat our calculations at higher resolu- 
tion and with more realistic initial gas density and temper- 
ature profiles (adapted from a standard a-disk). In § |3.5| we 
discuss possible numerical issues, and in § |3.6| we compare 
our results to previous work. Finally, in § El we summarize 
our main conclusions and their implications, and outline nat- 
ural future extensions of this study. 



© 0000 RAS, MNRAS 000, 000-000 



Gas Disk Response to Black Hole Recoil and Mass Loss 3 



2 INITIAL CONDITIONS AND SIMULATION 
DETAILS 

Our calculations were performed with version 2.5 of the pub- 
licly available Advanced Simulation and Computing (ASC) 
Flash code, developed at the University of Chicago. This 
code uses the split piece- wise parabolic method (PPM), 
making it well-suited for modeling shocks ( jFryxell et al.| 
20001. We use a uniform-mesh grid in polar coordinates to 
simulate a two-dimensional vertically averaged disk. Our 
runs cover a physical duration much shorter than the vis- 
cous time-scale, and we therefore do not include any phys- 
ical viscosity. The central black hole is chosen have a mass 
of M = 1O 6 M0, which is approximately the most favorable 
value for detection by LISA (although our results can be re- 
scaled to apply to BHs with different masses; see § 3.5 1. The 



inner edge of the grid is located 100R a (Schwarzschild radii), 
the radius at which the central gap opened by a binary black 



3.5 



hole is expected to freeze ( jMilosavljevic fc Phinney 2005 
|Haiman et al.|2009b[ ). 

A high-resolution mesh is used for radii between 100 < 
R/R a < 10 4 for all simulations, except the runs Kadl-5 (see 
Table [l] below), which use a high-resolution mesh between 
100 < R/R s < 1.5 x 10 4 . A lower-resolution mesh extends 
to 3.1 x 10 5 i? s . We tested the effect of boundary conditions 
by simulating a disk in Keplerian motion without the influ- 
ence of mass loss or kick. Within the 486 days we simulate 
after the merger, numerical fluctuations contributing to the 
density and luminosity of the disk were confined within the 
100 < R/R s < 10 3 region. The amplitude of the numeri- 
cal density perturbations was well below the over-densities 
produced by shocks in the disturbed disks, so in our anal- 
ysis of the densities and morphologies of the shocks, we in- 
clude the entire high-resolution mesh region, down to 100R a . 
However, the contribution of this low-amplitude noise to 
the luminosity was still significant. As a consequence, our 
post-processing analysis related to the disk luminosity is re- 
stricted to radii between 10 3 < R/R s < 10 4 ; corresponding 
to 30 days < t < lyear (see discussion of this issue in 
below) . 

In our fiducial run (see "Kiso3" in Table [I] below) , the 
high-resolution grid contains zones with sizes Ar = 20i? s 
and A9 = 0.01 radians, which is the effective spatial reso- 
lution of this run. This run had a corresponding numerical 
time-resolution, in physical units, of At « 35 seconds. How- 
ever, to save disk space, FLASH outputs were made only 
at every w 10 4 th time-step, allowing us to sample the re- 
sults on time-scales of ~ 4 days. In physical units, the total 
simulated time interval is 486 days (or ~ 1.4 X 10 7 M, in 
gravitational units with G = c = 1). For comparison, the 
orbital time at the inner and outer edge of our disk is ~ 4 
hours and 0.5 years, respectively, so that our simulations 
cover approximately 3000 and 3 full orbits at these radii. 
To test the numerical convergence of our results, we per- 
formed three additional variants of the fiducial run, one at 
lower resolution, with (Ar, A9) — (40 R s , 0.02 rad), and two 
at higher resolutions of (Ar, AO) — (10R S , 0.005 rad), and 
(Ar,A8) = (5 R s , 0.0025 rad). The corresponding physical 
time-steps in these runs were At = 70, 17.5, and 8.75 sec- 
onds, respectively. 

It has been shown that the black hole binary can ex- 
cite eccentricities in the disk within a region comparable 



to the binary separation (MacFadycn fc Milosavljevic|2008 



Hayasaki et al.|2007||2008||Cuadra et al.|2009[ ). This region 
shrinks as the binary separation shrinks, and our region of 
interest (> 10 3 R S ) is an order of magnitude larger than the 
binary separation when gravitational radiation first domi- 
nates the inspiral timescale. We therefore assume that the 
gas orbits have had time to circularize, and initially follow 
zero-eccentricity Kepler orbits. 

The disk is assumed to consist of ideal gas, and for sim- 
plicity, we use the constant gamma-law equation of state 
employed by FLASH. We chose the value 7 = 5/3, corre- 
sponding to monatomic ideal gas. However, the equation of 
state can be calculated using one of two methods. The first 
method is used in simulations we call "adiabatic," which use 
the internal energy of each mesh zone to calculate the tem- 
perature and pressure for the respective zone. The second 
method is "isothermal" for which temperature is used to cal- 
culate the internal energy and pressure for each zone. The 
temperature used in these calculations is held fixed through- 
out the disk for the duration of the simulation at a pre- 
specified value. In practice, this forces P oc p, but keeps 
Ei nt constant. Our adiabatic and isothermal runs should be 
considered limiting cases, which bracket the true behavior 
of a thin disk with a more realistic treatment of radiative 
cooling. 

For simplicity, in a first set of runs, we assume the 
disk initially has a constant surface density (E = 1.5 x 
10 5 g cm -2 ; the total mass of the disk is 2.1 x 10 3 M Q ) and a 
constant temperature (chosen from the range 4.5 x 10 2 K ^ 
T ^ 1.4 x 10 s K). These choices are broadly motivated by 
the physical values expe cted in a standard a-disk around a 
M black hole (e.g. |Frank et al.||2002[ ). 



10 



In a second set 



of runs, we adopt radial power-law profiles for the surface 
density, temperature, and scale-height, given by 
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The scale-height H(R) is not actually used in the 2D simula- 
tion, but we use it in § |3.6| below to compute the 3D gas den- 
sity, needed to predict the Bremsstrahlung luminosity. For 
a standard a-disk, with common choices for the viscosity 
parameter, disk opacity, and other parameters, the power- 
law slopes of the physical parameters are somewhat different 
in the "inner" (radiation pressure- and electron scattering 
opacity-dominated) and "middle" (gas pressure- and elec- 
tron scattering opacity-dominated) regions. The transition 
between these two regimes occurs around ~ 10 3 - 10 4 i? s 
for typical disk parameters, but the power-law slopes in 
these two regimes are, in any case, quite similar (see, e.g., 
Haiman et al.||2009bl. The power-law slopes we adopt in 



~aLp 



equations |1|-([3J are averages between the values in these 
two disk regions. Note that the total mass of our a-disk is 
580 M , somewhat lower than the mass of our uniform disk. 

The presence of the circumbinary disk is expected to 
align the spins of both BHs during the early stages of their 
evolution (at large separations), to be parallel with the an- 



© 0000 RAS, MNRAS 000, 000-000 



4 Corrales et al. 



Table 1. Summary of constant initial surface density disk sim- 
ulation runs. In all runs, the central BH has a mass of 10 6 Mq , 
and the disk extends from 10 2 to 10 4 Schwarzschild radii. In runs 
that include a kick, the kick velocity is 530 km s — 1 , oriented in 
the plane of the disk, and in runs that include mass loss, the 
fractional loss is assumed to be 5%. 



simulations. The runs were performed at the Center for Cos- 
mology and Particle Physics at New York University using 
the Ria cluster. Ria has 368 2.6 GHz CPUs with Infinband 
interconnect and 768 GB of memory. For reference, we note 
that our fiducial run, labeled Kiso3, was run on three nodes 
(24 processors), and took four days to complete. 



Simulation Type 


Adiabatic 


Isothermal 


Temperature 


kick only 


Kadi 


Kisol 


450 K 


(no mass loss) 


Kad2 


Kiso2 


4500 K 




Kad3 


Kiso3 


1.4 x 10 5 K 




Kad4 


Kiso4 


4.5 x 10 5 K 




Kad5 


Kiso5 


1.4 x 10 6 K 






Kiso6 


1.4 x 10 s K 


Higher— resolution 








version of Kiso3 




Kiso3+ 


1.4 x 10 5 K 


Highest-resolution 








version of Kiso3 




Kiso3++ 


1.4 x 10 5 K 


Low-resolution 








version of Kiso3 




Kiso3~ 


1.4 x 10 5 K 



mass loss only 
(no kick) 
mass loss and kick 



Mad3 Miso3 1.4 x 10 5 K 

MKiso3 1.4 x 10 5 K 



Table 2. Summary of simulation runs with an initial ct-disk 
profile. As in the constant surface density runs, the central BH 
has a mass of 10 6 Mq, and the disk extends from 10 2 to 10 4 
Schwarzschild radii; the kick velocity is 530 km s — 1 , oriented in 
the plane of the disk, and the fractional mass loss is 5%. The res- 
olution in all a-disk runs is the same as in the highest-resolution 
constant-density run Kiso3 ++ . 





Simulation Type 


Abbreviation 


Isothermal 








kick only (no mass loss) 


Kiso 




mass loss only (no kick) 


Miso 




kick and mass loss 


MKiso 


Adiabatic 








mass loss only (no kick) 


Mad 




kick and mass loss 


MKad 



gular momentum of the disk (Bogdanovic et al. 2007[ ). In this 
case, the recoil will be oriented within the plane of the disk; 
in any case, the two-dimensional nature of our simulations 
precludes us from studying out-of-the pl ane kicks. We fur- 
ther use the typical value of 530 km s _1 ( Baker et al.|2008 1 
throughout this paper. Our simulations are set up in the ref- 
erence frame traveling with the kicked SMBHB. Since the to- 
tal disk mass within our simulated region of 100— 10, 000i? s 
is only a fraction ~ 10~ 3 of the SMBHB mass, we neglect 
the inertia of the gas bound to the SMBHB. In practice, 
we simply add a unidirectional velocity component to the 
gas, everywhere in the plane of the disk, in our initial condi- 
tions, on top of the circular Kepler velocities. In runs which 
include a mass loss, we adopt a fractional loss of 5% - i.e, 
the disk is set up as before, but the central mass is reduced 
to 9.5 x 10 Mq at the first time-step of the simulation. 
Tables [I] and [2] summarize the parameters of all of our 



3 RESULTS AND DISCUSSION 

3.1 Post-Merger Disk Evolution and Shock 
Structure 

Figure [l] shows snapshots of the two-dimensional surface 
density distribution in our fiducial isothermal (Kiso3; left 
panel) and adiabatic (Kad3; right panel) simulations t = 210 
days after the merger. As the figure shows, sharp overdensi- 
ties develop in both cases, with a morphology similar to the 
spiral caustics in the collisionless case ( |Lippai et aLl|2008[ ). 
The density contrasts are visibly higher in the isothermal 
case, as one intuitively expects from the higher compress- 
ibility of isothermal gas. 

The sharpness of the thin arcs traced out by the over- 
dense regions suggest the presence of shocks. To assess 
whether shocks have indeed formed, once the simulation 
runs have completed, the output files were tested with an 
algorithm identical to the shock_detect feature available in 
FLASH. To receive a shock flag, a zone must pass both of 
the following criteria. First, the pressure difference between 
the zone and at least one of its neighbors must exceed the 
difference expected from the Rankine-Hugoniot jump con- 
ditions for a shock of a pre-specified minimum Mach num- 
ber M. To allow for oblique shocks, while comparing each 
zone to its neighbors in different directions, the pressure is 
weighted according the fraction of the total velocity made 
up by the velocity vector in the appropriate direction. Sec- 
ond, to avoid falsely flagging expanding regions as shocks, 
the velocity divergence at the zone in question must be neg- 
ative. 

The zones identified by the above criteria shows that 
shocks indeed arise in our runs at early times, and propa- 
gate outward. Shocks at early times in the disk were found to 
be typically weak. For example, in our Kiso3 run, at t — 30 
days, among the shocked zones identified with the threshold 
M > 1.1, we find that 20% exceeded M > 2. Subsequently, 
the shocks increased in strength, as they moved outward 
(with the top 20% of the shocked zones having M > 3 by 
t — 200 days). Shocks appeared in all of our runs except 
Kiso6 - in this run, the constant isothermal temperature was 
increased to correspond to a sound speed of c a ~ 10 3 km s _1 , 
approximately two times larger than i>kick- This suggests 
that i>kick > c s can be used as an approximate criterion 
for the kick to cause shocks. This is consistent with the ex- 
pectation in collisionless disk, in which particles cross their 
epicyclic orbits with a relative velocity of Av ~ fkick/iWbit; 
in bound regions of the disk, this relative velocity is limited 



to Av < Wkick (Lippai et al. 12008 1. 

In Figure [2] we sfiow the radial position of the densest 
point identified at several different times in the isothermal 
(Kisol-6; upper panel) and adiabatic (Kadl-5; lower panel) 
runs, together with the surface density enhancement at this 
point, defined as AE = E/E; (where E< is the initial con- 
stant surface density). The points are shown every 4 days, 
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Figure 1. Contour images of the surface density distributions in our fiducial isothermal (Kiso3; left) and adiabatic (Kad3; right) 
simulations at t = 210 days following the merger and the recoil of the SMBHB. In both cases, the rotation of the disk is clockwise, and 
the direction of the kick is up along the y axis; the disks are shown in Cartesian coordinates extending to 10 4 i? s . The color scales at the 
bottom are in units of the overdensity, relative to the constant initial surface density. 



between t = 4 and t = 300 days after the merger, with 
4 days corresponding to the left-most point, and 300 days 
to the right -most point. In the isothermal runs the densest 
point in the shock wave tends to stay at the same azimuthal 
position in the disk, ~ 37r/2 rad (as measured clockwise 
from the y axis), which corresponds to the region of the 
disk where the kick direction is parallel to the tangential 
velocity of gas in the disk. The kick adds constructively to 
the Keplerian velocity of gas at the angle 3n/2, and it is not 
surprising that gas at this azimuthal position achieves the 
largest over-densities first. 

In Figure [3] we show the radial position of the densest 
point in the disk as a function of time, in the same 10 runs as 
in Figure [2] (Kiso6 is excluded). For comparison, the black 
squares show the location of the outermost caustic identified 
in a collisionless case ( |Lippai et al.|2008| see further discus- 



3.6 below). The location of the isothermal shocks 



follow the caustics quite closely. However, interestingly, af- 
ter t « 100 days, the densest shock appears to propagate ~ 
twice as rapidly in the adiabatic simulations. As Figure [I] 
shows, the overall size and morphology of the sharp spiral 
patterns forming in the isothermal and adiabatic runs are 
quite similar, except that in the adiabatic runs, the spiral 
arms are more diffuse and each arm tends to extend to a 
somewhat larger radius than in the corresponding isother- 
mal case. The difference in the apparent radial propagation 
speed of the densest point arises because typically, in the 
adiabatic runs (i) the densest point falls further out along 
the outermost spiral arm, and (ii) the outermost spiral arm 
extends to slightly a larger radius. Furthermore, the dens- 
est point is typically farther out along the spiral pattern in 
adiabatic disks with increasing temperatures. Interestingly, 
the adiabatic run with the hottest temperature (Kad5) is an 
exception; in this run, the densest point again falls further 



in along the outermost spiral arm, which is noticeably more 
diffuse compared to the Kad4 run. The spiral arms in the 
adiabatic runs also tended to diffuse over time, with pairs 
of adjacent arms eventually overlapping. In comparison, the 
isothermal spiral arms remained distinct and did not inter- 
act. This likely explains why the isothermal curves are much 
smoother than the adiabatic curves. 

These figures also illustrate the strength of the den- 
sity perturbations caused by the kick, allowing us to draw 
several interesting conclusions. In particular, as the tem- 
perature, and therefore the pressure increases, the density 
enhancements generally decrease, as expected. The two low- 
est temperature cases (T = 450 and 4500K) differ very little 
in either the adiabatic or isothermal case, indicating that 
pressure effects are negligible when the disk temperature is 
below T « 4500-K\ Also as expected, the overdensities are 
much larger in the isothermal runs, reaching AE > 40 in 
our fiducial Kiso3 case, and still rising at the end of the 
run, whereas in the adiabatic run with the same tempera- 
ture (Kad3), the overdensity reaches a plateau and is limited 
to AE < 10. Pressure has little effect at very early times 
(except for the extremely hot disk Kiso6), when the shocks 
are within < 10 3 i? s and have relatively low-overdensities. 
Note that the sound crossing time near the inner edge of 
the disk in our fiducial run is R/c s ~ 70 days, comparable 
to the time when pressure starts to have an impact on the 
shock overdensities. Interestingly, pressure suppresses den- 
sity enhancements at late times by large factors. This last 
conclusion is somewhat surprising, given expectations from 



3.6 



collisionless disks (see further discussion of this point in 
below) . 

In Figure|4] we show the total mass of the shocked mate- 
rial divided by the total mass of the disk at a given time, all 
within the high-resolution mesh region (100 < R/R s < 10 4 
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Figure 2. Evolution of the maximum surface density in our 
isothermal (Kisol-6; upper panel) and adiabatic (Kadl-5; lower 
panel) runs. The points show the radial position and the local 
overdensity at the densest point in the disk, identified at several 
different times between t = 4 and t = 300 days after the merger 
(with 4 days corresponding to the left-most point, and 300 days 
to the right-most point). 
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Figure 4. Mass fraction of the disk material undergoing a shock 
with a minimum Mach number of Ai 1.1, as a function of time. 
The different symbols correspond to the runs shown in Figure [3] 
The curves begin to flatten when the outermost part of the spiral 
pattern begins to leave the high-resolution mesh region, in both 
the isothermal (~ 250 days) and the adiabatic simulations (~ 400 
days; note that in the adiabatic runs, a larger mesh was used, 
extending to 15.000.Rs). 



for Kisol-5, and 100 < R/R a < 1.5 x 10 4 for Kadl-5) us- 
ing the shock criteria denned above, and a minimum Mach 
number of M — 1.1. For cold isothermal disks, a significant 
fraction of the disk (~ 20%) can be experiencing a shock at 
a given moment. The shocks in the adiabatic cases comprise 
a smaller fraction of the disk, because the temperature and 
pressure (oc pT) both increase in the compressed regions, 
reducing the bulk velocities and shock strengths. Overall, 
compared to the isothermal cases, in the adiabatic runs the 
density enhancements are reduced, and shocks are weakened 
and are less prevalent - these trends are as expected, and 
follow from the lower compressibility of adiabatic gas. 
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Figure 3. The figure shows the radial location of the densest 
point in the disk as a function of time. The same runs are shown 
as in Figure[2] with the same notation (except Kiso6 is excluded, 
since it had no shocks). For comparison, the black squares show 
the location of the outermost caustic identified in a collisionless 
disk i ]Lippai et al,|2008^ . The straight dotted lines correspond to 
propagation at constant speeds of i;idck = 530 km s -1 (lower line) 
and 2^i t ; c ] I = 1060 km s _1 (upper line). 



3.2 Impact of Mass-Loss 

The results in the previous section included only the ef- 
fect of the kicks, and not the effects of the accompany- 
ing mass loss. Upon coalescence of the SMBH binary, en- 
ergy is carried away by gravitational waves, causing an in- 
stantaneous reduction of the rest mass by the amount of 
(Mi + M2)£gw- The fractional loss depends on the mass ra- 
tio, on the orbital parameters, and on the magnitude and 
orientation of both BH spins, for which the full param- 
eter space has not yet been explored in direct relativis- 
tic numerical simulations. For nearly equal mass binaries 
(q = M1/M2 ~ 1), however, the mass loss can be sev- 



eral percent (e.g. Tichy & Marronctti 2008 Herrman et al 
|2007a[ ). For simplicity, throughout this paper, we adopt a 
constant fractional mass loss of 5%. We initialize the disk 
with Keplerian velocities for the original pre-merger point 
mass (Mi + M2) = 10 6 Mq. However, in the remainder of 
the simulation, the gravitational potential is that of a central 
point mass with (Mi + M 2 )(l - e G w) = 9.5 x 10 5 Mq. 

The effects of such an instant mass-loss was treated re- 



cently by O'Neill et al. (20091 and Megevand et al. (20091 



in three-dimensional adiabatic disks. Here we extend these 
works by following the impact of the mass-loss case in an 
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Figure 5. Top panel: The left-hand side contour map shows the surface density for an isothermal 1.4 X 10 5 K disk, where the central 
BH experiences sudden mass— loss by 5% (Miso3). No kick is included. Concentric rings of density spikes develop and propagate outward, 
here shown at t = 210 days after the merger. The figure on the right shows the radial profile of the disk surface density for the isothermal 
disk (Miso3) in comparison to surface density from the adiabatic disk (Mad3) which was subjected to the same initial setup. Bottom 
panel: The contour map on the left shows the surface density of a disk at t = 210 days after the merger. A kick, with = 530 km s _1 

is included in addition to the mass loss (MKiso3). In the right panel, we also show, for comparison, the radial overdensity profile in 
our isothermal disk that includes only a kick, and no mass loss (Kiso3). Interestingly, the addition of the mass-loss results in a small 
reduction of the kick-induced overdensities (as well as a slightly faster outward-propagation of the outermost shock). 



adiabatic (Mad3) as well as in an isothermal disk (Miso3). 
In the latter case, we also investigate the combined impact 
of mass loss and recoil (MKiso3). As noted in the Introduc- 
tion, having one less dimension allows our simulations to 
cover the disk out to a much larger radius (10 2 — 10 4 R S 



O'Neill et aL|2009fr, and for a m uch 



5 hrs in |Q'Neill et al.|2009| for 



the range 2 — 50R 
longer duration (~ lyr vs 
a 10 6 M BH). 

The top portion of Figure[5]shows a surface density map 
of an isothermal disk, in the mass-loss case, around t — 210 
days after the merger (in the left panel). In the right panel 
of the same figure, we show the radial surface density profile 
of the same disk, together with the profile of an adiabatic 
disk with the same mass-loss. When the impact of the kick 
is not included, the perturbations are azimuthally symmet- 
ric, giving rise to concentric, outward-propagating ripples. 
The basic mechanism behind creating the ripples can be vi- 
sualized simply as follows: at the moment of mass loss, each 
particle, orbiting at its original Keplerian velocity, is moving 
at a speed that is too fast for a circular orbit around the re- 



duced central mass. As a result, all particles find themselves 
at the pericenter of a new, elliptical orbit. Summing the ef- 
fective contributions of each particle to the surface density, 
spread over these orbits, results in a net dilution of the disk. 
Furthermore, the radial epicyclic oscillations of the set of 
particles at a fixed initial radius will occur at a frequency 
that monotonically decreases radially outward, resulting in 
a pattern of quasi-periodic compression / decompression in 
the radial direction. The resulting density spikes have lower 
amplitudes, but are more frequent toward the inner bound- 
ary. These effects are clearly seen in Figure [5] except that, 
notably, the pressure forces in the adiabatic case are suffi- 
cient to almost completely wash out the innermost density 
ripples, out to R < 2, 000 R B , by t = 210 days. 

When a kick is included in addition to the mass-loss 
(MKiso3), the density enhancements are often lower in com- 
parison to the kick-only case (Kiso3). The bottom portion 
of Figure [5] shows a map of the surface density and a radial 
slice of the surface density profile along the azimuth 8 — n 
(running in the direction of the -y axis) in both cases. The 
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density enhancement from the mass-loss + kick run is ap- 
preciably lower than in the kick-only case, especially for the 
outermost, densest spike, at R m 4, 500 — 5, 000-R s . 

3.3 Observational Signatures 

Our adiabatic runs conserve energy, so, by definition, no 
energy loss occurs in these runs, and the disk has strictly 
zero luminosity. In reality, the disk will radiate - indeed, 
a prediction of its time-variable light-curve, following the 
merger, will be invaluable for identifying EM counterparts 
to LISA sources. A realistic prediction requires solving for 
the vertical structure of the disk, as well as a treatment 
of radiative cooling and radiative transfer across the disk 
(which is likely to be optically thick; see below). While this is 
not possible in our 2D simulations, our isothermal runs allow 
us to calculate the energy dissipation occurring in the disk - 
this can serve as a rough, order-of-magnitude estimate for 
the disk luminosity. 

More specifically, at each time-step, FLASH solves Eu- 
ler's equations, expressed in conservation form, to update 
the total specific energy (E) of each zone which includes the 
combined kinetic and internal energy. In our case, the disk is 
cold, and the internal energy (E- mi ) is very small compared 
to the kinetic energy of orbital motion, and is therefore a 
small fraction of the total energy E. This can cause numeri- 
cal inaccuracies, and we found that it was necessary to solve 
for _Bi n t explicitly, using an internal energy advection equa- 
tion that is separate from the total energjj^] The pressure 
(and temperature, if required) of the zone is subsequently 
obtained using the specified equation of state. 

When an isothermal equation of state is used, the pre- 
specified constant initial temperature T lso is used to cal- 
culate the internal energy, so that it is reset to E{^ = 
(7 — 1) _1 pkT m ° I pm p cc p. The difference of the internal 
energies before and after this re-set, AE- m t — -Bint — E^, cor- 
responding to the pAV work done by the gas, is effectively 
radiated instantaneously; we compute the corresponding lu- 
minosity by dividing A_Bi nt by the simulation time-step. 

The light curve obtained by this method corresponds to 
the limiting case of efficient radiative cooling that keeps the 
gas at constant temperature. An ambiguity in this interpre- 
tation arises, however, in zones which expand, and therefore 
would cool, in the adiabatic case. The isothermal condition 
requires that these fluid elements gain, rather than loose en- 
ergy, in order to maintain their constant temperature. The 
above calculation simply returns a negative luminosity in 
these zones, which should either be excluded from the total 
luminosity budget (if, in reality, these zones remain cooler 
than T lso ), or may be included (under the assumption that 
the corresponding heating comes at the expense of the ra- 
diation generated in the compressed regions). Clearly, the 
isothermal approach is not adequate by itself to resolve this 
ambiguity, and below we will simply quote luminosities both 
with and without the negative contribution of these expand- 
ing zones. These may also be regarded as upper and lower 
limits; our light-curves are inherently limited to an accuracy 

1 See equation 13.7 in the FLASH User 
Guide, Version 3.2, July 2009; available at 
http:/ /flash. uchicago.edu/website/codesupport/flash3 _ug_3p2.pdf 
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Figure 6. A snapshot of the luminosity density, defined as the 
the energy radiated per unit area and unit time, in our fiducial 
Kiso3 run, here shown at t = 210 days. The bright regions, with 
a high positive luminosity, closely trace the spiral shocks seen in 
Figure [T] The color scale shown at the bottom of the figure is in 
units of erg s" 1 cm -2 . Note the negative contributions from the 
post-shock regions (see text for discussion). 



comparable to the difference between these two estimates 
(which we find is typically a factor of ~ two) . 

Figure [6] shows a snapshot of the luminosity density, de- 
fined as the energy radiated per unit area in a single time- 
step At, j = AEi n t/AAAt in our fiducial Kiso3 run, at 
t — 210 days. As noted above, the luminosity in the region 
inside R < 10 3 i? s is dominated by numerical noise. To be 
conservative, we have therefore excised these inner regions 
from our luminosity calculations. Outside R < 10 3 R S , the 
bright regions, with a high positive luminosity, clearly trace 
the shocks quite well, which is indeed what one expects, 
since these are the locations where the gas is being most 
compressed. Adjacent to the bright shocks, the figure also 
shows regions of post-shock decompression, where the ve- 
locity divergence is positive, and the luminosity is negative. 

In Figure [7j we show the corresponding light curves in 
our fiducial run (Kiso3) - again, conservatively, limited to 
the region 10 3 < R/R s < 10 4 - when only the positive- 
luminosity zones are considered (black solid curve) and when 
the contribution from the negative-luminosity zones is sub- 
tracted (black dotted curve) . The other curves in this figure 
show how the light-curves depend on numerical resolution, 
and will be discussed in § |3.5| below. The absolute value 
of the luminosity is significant, corresponding to w 16% of 
the Eddington luminosity LEdd for a 1O 6 M0 black hole in 
the first case, and reduced by less than a factor of two, to 
~ 10% I/Bdd in the latter case. We note that despite the 
significant luminosity, the total energy radiated in the first 
year after the merger (integrating the light curve in Fig- 
ure M) is found to be ~ 0.1 x Afdisk^kick- In principle, the 
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Figure 7. Light-curves obtained in T = 1.4 X 10 5 K isother- 
mal disks experiencing a kick, by integrating the luminosity den- 
sity, such as shown in Fig. [6] across the surface of the disk. The 
solid curves include only the positive— luminosity zones, whereas 
the dotted curves include the contribution from the negative- 
luminosity zones. The black curves are derived from our fiducial 
run (Kiso3), and the magenta/rcd/blue curves show the results 
of a resolution study (corresponding to the high-resolution runs 
Kiso3 ++ and Kiso3 + ; and the low— resolution run Kiso3~ , respec- 
tively). The total positive luminosity in the highest-resolution run 
(top curve) reaches Ri20% of the Eddington limit for a 10 6 Mq 
black hole. 
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Figure 8. Contributions to the total (positive) luminosity from 
different annuli in the fiducial run Kiso3, as labeled. The densest 
shocked region of the disk tends to dominate the total luminosity, 
and each annulus provides a maximum contribution when this 
"bright spot" moves across it. 



kick-induced dissipation could extract and radiate away a 
non-negligible fraction of the total potential energy of the 
disk; the luminosity we find implies an energy release well 
below this upper limit. 

In Figure |SJ we decompose our light-curve, and show 
the contribution of the positive-luminosity zones from differ- 
ent annuli in the disk, binned by radius. The densest shocked 
region of the disk tends to dominate the total luminosity, 
and, as the figure shows, each annulus provides a maximum 
contribution when the "bright spot" moves across it. The 
luminosity rises steadily overall, which is attributable pre- 
dominantly to having more radiating material at larger radii 
(the bins shown in Fig. [8] have equal radial widths, so their 



•p 1.0x10"- 




300 
erqer (days) 



Figure 9. Light curves as in Figure [7] except here we contrast 
the results from the mass- loss only (Miso3), mass-loss + kick 
(MKiso3), and the kick-only cases (Kiso3). When the effect of 
mass-loss is included in addition to the kick, it reduces the overall 
luminosity of the disk by a modest factor. In the case of mass- 
loss only, the luminosity of the disk is reduced significantly (by a 
factor of fa 4). 



area grows linearly with radius). As discussed above, the 
shocks also get stronger with time; however, this is a rela- 
tively weak effect, and the overall rise in the luminosity, seen 
in Figure]?] roughly tracks the build-up of the shocked mass 
(see Fig. El . At the most luminous portion of the shock, in 
our fiducial run, the temperature typically rises during each 
time step by « 4%, before it is re- set to 1.4 x 10 s K. 

When the effects of mass-loss are included in the sim- 
ulation, in addition to the kick, the luminosity of the disk 
tends to decrease. This can be expected from Figure[5] which 
showed that the overdensities in the spiral shock waves tend 
to be reduced compared to the case without mass-loss. Fig- 
ure [9] illustrates the effects of mass-loss on the light curve 
produced from a kicked disk, as well as the luminosity pro- 
duced in the mass-loss case. 

The above estimates for the light-curve utilize only the 
effective dissipation that occurs, implicitly, in our simula- 
tions. The photons generated by this dissipation must escape 
the disk before they can be observed. We next use standard 
a-disk models to estimate the basic parameters of a more 
realistic disk. We adopt the equations describing the radial 
profile of the density, scale-height, temperature, and opacity 
of a standard thin accretion disk as summarized in IHaimanl 
et al. ( 2009b based, in turn, primarily on Frank et al. 2002 



and on Goodman & Tan 2004). In Figure 10 we show the 
temperature (T; in Kelvin; green curves), the vertical op- 
tical depth (r; blue/red when electron scattering/free-free 
opacity dominates, respectively), the cooling time due to 
Bremsstrahlung emission (t coo i; taken from eq. 6.53 in |Frank] 
|et al.|2002| in units of seconds; magenta curves) and the ver- 
tical photon escape time (t esc ; in units of days; black curves). 
The vertical optical depth is taken to be r = max(r cs , ), 
where r cs is the free-free optical depth from Kramer's opac- 
ity, and t os is the electron scattering optical depth. The 
photon escape time is obtained as the product of the mean- 
free-path H/t (where H is the vertical scale-height) and 
the number of scatterings in a random walk required to es- 
cape (r 2 ), divided by the speed of light, t CBC = (H/c)t. All 
curves are for a M = 1O 6 M0 black hole. Solid curves are 
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Figure 10. Physical quantities in a standard a-disk around 
a M = 10 6 Mq black hole, as a function of radius (lower la- 
bels) and corresponding orbital velocity (upper labels). We show 
the temperature (T; in Kelvin; green curves), the vertical opti- 
cal depth (t; blue/red when electron scattering/free— free opacity 
dominates, respectively), the Bremsstrahlung cooling time (t coo i; 
in units of seconds; magenta curves) and the vertical photon es- 
cape time (t CBC ; in units of days; black curves). Solid curves are for 
a standard a-disk (with viscosity scaling with total gas + radia- 
tion pressure), and the dashed curves are for a standard /3-disk 
(with viscosity scaling with gas pressure. The thin vertical lines 
bracket the extent of the high-resolution mesh in our simulations. 



for a standard a-disk (with viscosity proportional to the to- 
tal pressure of gas + radiation), and the dashed curves are 
for a /3-disk (with viscosity proportional to the gas pres- 
sure). The viscosity parameter, the accretion rate (in Ed- 
dington units), and the radiative efficiency were chosen to 
be a = 0.3, m = 0.1A#Edd and e = 0.1, respectively. All 
the other parameters, which have a relatively minor effect 
on the profiles, are the same as the fiducial ones defined in 



Haiman et al. (12009b I 



As Figure |10| shows, the disk is optically thick (r ^> 1) 
at all relevant radii, and at radii < 80007? s , the opacity is 
dominated by electron scattering. Near this radius, the lo- 
cal Bremsstrahlung cooling time is negligibly short, and the 
photon escape time from the disk is only ~ 3 days. There- 
fore, even though the disk is optically thick, the photons, 
generated by dissipation in the mid-plane, will emerge from 
the disk with a relatively short delay. These conclusions 
are derived for a steady disk. In the case of a binary, the 
pile-up of material near 100R S will modify the disk struc- 
ture, corresponding, locally, to a non-accreting (M = 0) 
disk ( Ivanov, Papaloizou, fc Polnarev[T9 99 Milosavljevic & 



Phinney 2005 I; however, the steady-disk model should give a 
good order-of magnitude estimate farther out at ~ 10007? s , 
where the disk is relatively less disturbed before the merger. 

Since the disk is very optically thick, the spectrum 
that emerges will be significantly processed, and close to 
a black-body shape (sometimes referred to as "grey-body" ; 
e.g. Milosavljevic fc Phinney|20 05). Assuming for simplicity 
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Figure 11. Top panel: effective black-body temperature T c g for 
the most luminous point in the disk in the fiducial Kiso3 run, 
defined by requiring the black-body surface brightness crT^j to 
equal the dissipation rate derived in the simulation. Bottom panel: 
total monochromatic luminosity emerging from our fiducial Kiso3 
disk at three different times, obtained by integrating the local 
black body emission over the face of the disk. The three curves 
(from bottom to top) show the spectra at t = 80, 240, and 400 
days after the merger. 



that the radiation escaping from the surface has a perfect 
thermal spectrum, we can compute the required effective 
temperature T e g, by equating the black-body luminosity 
per surface area aT^ s to the luminosity density j we de- 
rive from the simulations. In the top panel of Figure |11[ 
we show the evolution of this effective temperature, com- 
puted as a function of time at the most luminous point in 
the disk. As the figure shows, T e g rises and asymptotes to 
~ 4.5 x 10 4 K, a value that is a factor of ~ 3 lower than the 
constant isothermal temperature T ISO imposed on the disk. 
This is consistent, to within a factor of ~two, with the ver- 
tical structure in thin accretion disk models, in which the 
surface temperature is lower by a factor of ~ r -1 / 4 than the 
mid-plane temperature (e.g. Frank et al. 2002). 

From the similarity of the effective temperatures, we 
infer that the transient dissipation required to maintain a 
constant temperature in the disturbed disks is comparable 
(albeit an order of magnitude higher) than the viscous dis- 
sipation that occurs in a standard, steady-state, thin accre- 
tion disk. This implies that if the pre-merger disk luminosity 
was comparable to that of a thin accretion disk, the merger 
would be accompanied by an order-of-magnitude rise in lu- 
minosity. In practice, the total luminosity of a steady disk 
is dominated by its inner regions, and the pre-merger lu- 
minosity of the circumbinary disk can be significantly lower 
(and correspondingly, the merger will result in brightening 
by a much larger factor). 

Interestingly, the monotonic rise in T c ff suggests that 
the shock-induced emission may have the unique signature 
of hardening in frequency by a factor of ~ two (in addi- 
tion to monotonically rising in luminosity by a factor of a 
few) within the year following the merger event. In the bot- 
tom panel of Figure |11| we show the emergent composite 
spectrum, integrating the local black body emission over the 
face of the disk, which shows the same two features: the disk 
brightens, and the spectrum hardens with time. If confirmed 
by more accurate modeling of the emergent spectrum, this 
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Figure 12. Total luminosity (black dotted line) compared to the 
luminosity contribution from different annuli (solid colored lines) 
in an a-disk that received a 530 km s _1 kick and experienced 5% 
mass loss from the central black hole (MKiso). The initial flares 
occur when the first two arms of the shock pass through the disk. 
However, the luminosity of the outer region of the shock decreases 
dramatically because less gas is swept up by the shocks in the 
outer regions of the disk, which is more diffuse. A large number 
of thin, tightly wound density waves persist throughout the entire 
disk, and the disk maintains a nearly constant luminosity until 
the end of the simulation. 
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Figure 13. Light curves for three isothermal disks, experiencing 
kick and no mass-loss (Kiso), mass-loss and no kick (Miso), and 
both kick and mass— loss (MKiso) as in Figure|9] except the disks 
here have initial power— law density and temperature profiles, ap- 
propriate for a standard «-disk, rather than constant values. The 
luminosities rise more rapidly, due to the higher initial density 
concentration near the central regions, but settle at values sim- 
ilar to those reached in the constant density and temperature 



should greatly help the identification of EM counterparts, 
since other types of known transient UV/optical/IR sources 
do not have these qualitative features. 



3.4 Varying the Initial Disk Profile 

All results in the previous sections are derived, for simplicity, 
for disks with constant initial density and temperature. In 
order to quantify the sensitivity of our results to the initial 
profiles, we have re-run several of our simulations with the 
power-law profiles in equations Q-|3j|, approximating those 
in a more realistic, standard a-disk. The spatial resolution 
of all the a-disk run is set to be the same as in the highest- 
resolution constant density case, Kiso3 ++ , with (Ar, A8) — 
(5 R s , 0.0025 rad). Table [2] summarizes the a-disk runs. 

Figure [12] shows the total luminosity of the a-disk, 
again conservatively from the restricted region 10 3 < 
R/R s < 10 4 (black dotted line) together with the luminosity 
contribution from different annuli (solid colored lines). The 
first peak in each annulus, which is easily identifiable around 
t = 10 days in the 1000 < R/R s < 2000 bin, comes about 
when the first arm of the spiral shock forms. The second 
winding on the spiral arm, which is considerably narrower 
at early times than the first spiral arm, accounts for the 
second rise in the binned luminosity (~ t — 50 days in the 
first bin). As the spiral pattern progresses outward, tighter 
windings of the spiral pattern persist and are generally long 
lasting due to the lack of viscous dissipation in the disk. As 
a result, after the first few windings of the spiral arm pass, 
each annulus contributes a steady background luminosity 
due to the numerous spiral density ripples that persist to 
the end of the simulation. While this is true in the uniform 
disk runs, as well, here the innermost annuli continue to 
dominate the total luminosity at all times, even though the 
outer annuli cover a larger surface area of the disk. It is 



also noticeable that, as the first few spiral arms propagate 
outward, their contribution contributes less and less to the 
total luminosity (compare Figure 



12 



to 



3 



where the flat 



portion of the lightcurve-contribution from each bin is due 
to the densest point moving into the respective region.) In 
the a-disk simulation, the front portion of the spiral shock 
contributes less to the total luminosity as it moves outward 
because the surface density of the gas is considerably lower 
at larger radii. 

Figure[l3]shows the light-curves in the three isothermal 
a-disk simulations. The absolute value of the luminosities 
arising from shocks in the a-disks is very close to the pre- 
vious results in the constant density disk cases. However, 
there are a few qualitative differences worth noting. First, 
the luminosity rises more rapidly in the a-disks, reaching 
peak brightness in m 100 days, compared to the constant- 
density disks, in which the peak is reached only after w 300 
days. While some of this more rapid rise can be attributed 
to the higher resolution employed for the a-disks, the res- 
olution study in Figure [7] makes it clear that most of the 
difference is caused by the change in the profiles. This qual- 
itative difference is unsurprising, since the a-disks are more 
centrally concentrated, and have more material in the in- 
ner regions, increasing the luminosities at early times. The 
mass-loss only case shows that a brighter and more prompt 
flare erupts after the merger, compared to the constant den- 
sity disk; again this is attributable to the larger amount of 
material near the central regions. When both mass-loss and 
a kick is included, we find a prompt flare, followed by a 
nearly steady luminosity, reduced by around 15% from the 
kick-only case. 

Overall, the light-curves in the a-disk and constant- 
density cases are reassuringly similar, indicating that vari- 
ations in the assumed initial disk profiles do not lead to 
order-of-magnitude changes in the predicted brightness of 
the after-glow. 
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3.5 Numerical Issues 

In this section, we discuss numerical effects from the inner 
and outer boundary of the simulation box, the dependence 
of our results on the spatial and time resolution, as well as 
numerically scaling our results to apply to black holes and 
disks with different physical sizes and kicks with different 
speeds. 

Boundary conditions. Both at the inner and the outer 
boundaries, we imposed "outflow" boundary conditions that 
maintain a zero pressure and density gradient at the edges 
of the simulation. To test the effect of these boundary con- 
ditions, we checked whether an unperturbed disk (without 
any kick or mass-loss) maintains its uniform circular Kepler 
motion. We found that the outer disk, beyond R s = 10 3 R S , 
indeed remained undisturbed, but modest density ripples 
(with over-densities E/Ej < 2) were produced artificially 
inside 10 3 R a . The amplitude of these numerical fluctuations 
are below the density enhancements excited by shocks in 
the disk (see for example Figures [5] and [3]) . Nevertheless, 
this low-amplitude noise covers a large fraction of the sur- 
face area of the inner disk, and can dominate the total lu- 
minosity in the inner regions. This is a concern especially 
for the q disks, in which the inner regions were also found 
to dominate the total luminosity. More specifically, we com- 
puted the luminosity inside R s — W 3 R S in the unperturbed 
constant density and a-disk models, and we found that they 
reach values larger than (in the constant surface density disk 
cases) or comparable (in the a-disk cases) to the luminosi- 
ties exhibited in the perturbed disks. In future work, we will 
attempt to identify the origin of this numerical noise, and 
to reduce its contribution to the luminosity. In the present 
paper, we simply derived luminosities outside R > 10 3 R S , 
which are safe from numerical noise, but represent an un- 
derestimate the true total luminosity. 

A separate concern is that due to the finite size of our 
simulation, the outermost regions of the spiral shock wave 
eventually reach the outer boundary of our region of interest. 
In the Kiso3 run, the outer edge of the high-resolution mesh 
is reached at t ~ 250 days. Although the spiral shocks are 
subsequently still present farther out in the low-resolution 
mesh, we have excluded these poorly resolved shocks from 
our analysis and results. For example, the dense spiral arms 
moving out of the simulation box produce the turn-over in 
the percentage of the disk undergoing a shock, shown in 
Figure |4| (the fi gure only includes shocked gas within 7? < 
10 4 R a ). 

Numerical resolution. In order to test the numerical 
convergence of our results, we re-computed the light-curves 
in our fiducial disk (Kiso3) at a factor of two lower (Kiso3~) 
and factors of two and four higher (Kiso3 + & Kiso3 ++ ) spa- 
tial resolutions. Note that FLASH chooses time steps adap- 
tively, matched to the spatial resolution, so each resolution 
test differs from the next by a factor of two in the time step, 
as well. The four curves in Figure [7] show the light-curves 
in each of these four runs, corresponding, from bottom to 
top, to increasing resolution. The light-curves begin to con- 
verge, in the sense that their fractional difference between 
pairs of resolutions begins to decrease. There is still a no- 
ticeable difference (by a relatively small fraction of « 10%) 
between the luminosities computed in our highest and next- 
to-highest resolution run. Given that the total luminosity 



is dominated by sharp features, it is not surprising that 
the disks are more luminous in the higher resolution runs 

- the density spikes are better resolved in these runs. At our 
fiducial resolution, the luminosities are underestimated by 
> 35%. Unfortunately, the computational time required to 
run the highest resolution simulations (which took approx- 
imately 2 weeks on 80 processors) was prohibitively expen- 
sive, and prevented us from further increasing our resolution, 
in order to achieve full convergence. However, the 10% dif- 
ference between the two highest resolution runs is modest 
in comparison to other uncertainties arising from the highly 
idealized nature of our luminosity estimates. The resolution 
study also shows that our somewhat under-resolved simula- 
tions are conservative, underestimating the true luminosity. 

Scaling simulation runs to physical units. The physical 
scales in each simulation are fixed by specifying the physi- 
cal values of any two independent dimensional parameters 

- these can be taken to be the BH mass M, and the kick 
velocity Ukick- Our runs can describe systems with different 
combinations of BH masses and kick velocities, provided the 
other parameters are suitably adjusted. Consider the quan- 
tity X, specified in the simulations by a dimensionless value 
A sim . The physical value A phys of this quantity can be ob- 
tained by a constant conversion factor /, and will scale with 
M and fkick as 

xph -=/(^r(^r. m 



10 6 M P 



V530 km/s/ ' 

Table [3] shows the factor fX slm for various parameters spec- 
ified in the initial condition of our fiducial run, and the 
mass- and velocity-dependence of these parameters; we also 
show the mass- and velocity-dependence of the physical lu- 
minosity and energy, derived from the simulations. For ex- 
ample, as «kick increases, « or b must increase in proportion, 
decreasing both the distance (oc v^? k ) and time (oc irf k ) 
scales. Thus, our results would be applicable, for example, 
for M = I0 6 Mq and v kick = I0 3 km s~\ but would then 
describe the evolution of a 4 times smaller disk, with an in- 
ner cavity of size 25R S , and for only « 50 days, rather than 
a year. The disks would also be 4 times hotter, and 16 times 
more massive than for the fiducial 530 km s _1 kick. The 
luminosity has the especially steep dependence oc v kick , so 
that higher-velocity kicks would be much more observable 
(provided, of course, that the scaled physical parameters are 
realized in nature). Other combinations to which our runs 
can correspond can be read off similarly from Table [3] 

3.6 Comparison with Previous Work 



Lippai et al. ([2008 ) considered the response of a collision- 
less disk to kicks. Qualitatively, the spiral shock structures 
shown in Figure [T] in our fiducial Kiso3 run are very similar 
to the caustics shown in Figure I of Lippai et al. ( 2008[ ) for 
a sy stem with the sam e parameters (the only difference is 
that Lippai et al.||200"8 used v k - lck = 500 km s _1 ). However, 
a more quantitative comparison is useful to separate hydro- 
dynamical effects from those arising from orbital dynamics. 

Figure [3] shows the radial position of the densest point 
in the simulated disks, compared to the location r c of the 



outermost caustic identified in a collisionless case (Lippai 
|et al.|[2008| ). The straight dotted lines correspond to prop- 
agation at constant speeds of u k ick = 530 km s _1 (lower 
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Table 3. Conversion from the physical parameters in our fiducial run to other values, for each parameter specified in our initial 
conditions, as well as for the luminosity and energy derived from the simulations. 



Simulation Parameter 


Conversion Factor (fX B1In ) 


7m 


Iv 


BH Mass (M Q ) 


10 6 


1 





Radius of inner disk edge (cm) 


3.0 x 10 13 


1 


-2 


Maximum time (days) 


500 


1 


-3 


Kick velocity (km s — *) 


530 





1 


Gas temperature (K) 


1.4 x 10 5 





2 


Surface Density (g cm -2 ) 


1.5 x 10 5 


-1 


4 


Luminosity (erg s _1 ) 







5 


Energy (erg) 




1 


2 



line) and 2«kick = 1060 km s _1 (upper line). As the fig- 
ure demonstrates, the shocks found in our isothermal runs 
form at locations very close to those of the collisionless caus- 
tics. Both follow the approximate simple linear expansion 



fkicki expected from the epicyclic approximation (Lip- 



pai et al. 2008 1, although the propagation is somewhat faster 



at early times, and somewhat slower at late times. Appar- 
ently, while the temperature has a strong impact on the 
overdensities and the strengths of the shocks (shown, e.g., 
in Fig. [5| it makes a comparatively small difference to their 
propagation speed or to their overall morphology. These re- 
sults also provide reassurance that the shock propagation 
was not affected by numerical edge effects from the inner 
and outer boundary of the simulation. Interestingly, on the 
other hand, the equation of state has a significant impact, 
with the propagation ~ twice as rapid in the adiabatic runs. 
As mentioned above, this apparent difference in the propa- 
gation speed of the densest point arises partly because the 
densest point falls further out along the outermost spiral arm 
in the adiabatic runs, and partly because the gas pressure 
makes the outermost spiral arms more diffuse and extend to 
a larger radius. 



As mentioned in § |3.1| above, we found that pressure 
suppresses density enhancements in the disk at late times, 
by large factors. This last conclusion is somewhat surpris- 
ing, given expectations from collisionless disks. Li ppai et al.| 
(2008) found that the relative velocity v c of particles, as 
they cross on their orbits to create the caustics, increases 
with increasing radius as v c ~ fkick/^orbit oc r 1 ^ 2 (where 
«orbit is the orbital velocity). The collisions therefore become 
super-sonic beyond some critical radius; this critical radius 
typically lies within the range of radii shown in Figure [2] 
For example, in the fiducial run Kiso3, the critical radius is 
at 1, 300-Rs, and in the higher-temperature runs Kiso4 and 
Kiso5, it is at 4, 2007? s and 13,000T? S , respectively. We find 
that the large overdensity of the shocks can be attributed 
to "pre-compression" . In particular, the shocks exceed the 
overdensity of a factor of 4 relative to the surface density at 
the start of the simulations. This is because significant over- 
densities start building up along the spiral pattern before 
the actual shocks occur. This increases the absolute density 
of the shock, since the shock is hitting gas that already has 
a surface density higher than the initial value. For example, 
in Figure [I] the adiabatic simulation (Kad3) shows that the 
outermost edge of the spiral pattern has overdensities form- 
ing with AE m 2. For a strong shock and 7 = 5/3, we 



expect a density enhancement about 4 times that of the ini- 
tial density, or E ~ 8Ei, which is what we see in the denser 
part of the spiral. A similar mechanism is at work in the 
isothermal simulation (Kiso3), but the shocked gas is able to 
cool and reach much larger densities. This pre-compression, 
which happens on time-scales much longer than the shock- 
formation itself, can apparently be significantly suppressed 
by gas pressure. 



Lippai et al. ( 2008 1 also proposed that the shock tem- 
perature scales with time/radius as T B h oc k oc v 2 oc v~^ it oc 
r c oc t c , and suggested that the luminosity is dominated by 
the outermost shocked shells, with both the total luminosity 
and the characteristic frequency of the emitted spectrum in- 
creasing with time. The same qualitative behavior is seen in 
our isothermal, constant-density simulations. However, we 
find that the largest Mach number in the spiral arms in- 
creases from M = 1.7 to M = 3, between t — 16 and 310 
days, much less steeply than linearly with time (and radius). 
In our adiabatic runs, the Mach numbers decline from 2 to 
1.1 in the same interval. This decline, however, is due to the 
increase in sound speed, rather than a decrease in the gas 
velocities. 

Assuming a disk profile with E oc r~ 3//5 (close to the one 
adopted in our a-disk runs), and also that a constant frac- 
tion of the disk mass interior to the outermost spiral caustic 
is shocked (M a hock oc Erdr) and that the energy is released 
over the time-scale of t c = r c /-Ukick, 
finds the scaling Lkick ~ (l/2)M s hock«c/tc oc r' 
cussed in § |3.3| above, the disk is likely to be optically thick, 
in which case the spectrum is significantly re-processed, in- 
validating naive conclusions about the optically-thin spec- 
trum. In the constant-density runs, we nevertheless find an 
approximately linearly rising luminosity, and also a mod- 
estly increasing characteristic frequency. While this is in 
qualitative agreement with the above expectations, most of 
the rise in the luminosity in our case is explained by the in- 
creasing fraction of the disk mass that is shocked, with the 
increase in the characteristic temperature contributing rel- 
atively little. In our a-disk runs, we find a light-curve that 
initially rises and then reaches a nearly constant luminosity 
after «100 days. 

As mentioned in the Introduction, |Q'Neill et al.| ( |2009| ) 
and Megevand et al. ( 2009 1 have recently employed three- 



Lippai et al. ( 2008 ) also 
kUc/t c oc r 1M/ ™~As dis- 



dimensional simulations to study the evolution of much 
smaller disks (extending from 2 — 50i? s and 10 — 30-Rs, re- 
spectively, compared to our 10 2 — 10 4 7? 3 ) on much shorter 
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time-scales (~ 5 hrs and ~ 6 hrs, scaling to a 10 6 Mq BH; 
compared to our ~ 500 days). Since our FLASH simulation 
results are output at times that sample the disk evolution 
only every « 4days, we can not directly compare the be- 



MKiso3 



havior of our disks to those in O'Neill et al. (20091 and 



Megevand et al. (20091 



Despite the mis-match in scales and different dimen- 
sionality, it is interesting to contrast the light-curves we 
find h e re wit h those in |Q'Neill et al] ( |2009| > and |Megevand| 



et al. (20091. We first note that the simulations in both 



O'Neill et at] (|2009|) and |Megevand et al.| (|2009[) are adi- 



abatic, i.e., they are energy-conserving by construction. 



O'Neill et al. ( 2009 I and Megevand et al. ( 2009 1 present esti- 



mates for the luminosity by computing the volume-integral 
of the quantity p 2 T x ^ 2 , which is proportional to the local 
Bremsstrahlung emissivity. As we discussed in § |3.3| above, 
the Bremsstrahlung cooling time is exceedingly short, so if 
dissipation indeed tracked the Bremsstrahlung emission, the 
adiabatic condition could not be maintained. If cooling were 
efficient, then one expects, instead, much stronger overden- 
sities to develop than in the adiabatic case, which would, in 
turn, generically increase the cooling rate tx j dV p 2 T 1 ^ 2 . 

Setting these inconsistency issues aside, in order to 
compare our results as directly as possible to O'Neill ct 
al. (20091 and Megevand et al. (20091, we computed the 



quantity I/Brom = J dVp z T 1/j: in both our adiabatic and 
isothermal simulations. Unfortunately, this comparison is 
further complicated by the fact that our simulations are two- 
dimensional. In order to calculate the 3D gas densities and 
volume elements, p and dV from the constant density disks, 
we assumed that the disk has the scale height H — c a /Q, 
obtained by assuming hydrostatic equilibrium in the ver- 
tical direction (here c a — \JkT '/ 'pm p and f2 = v^/r are 
measured directly in the simulation run). This assumption 
should be reasonable, since we simulated N 3> 1 dynami- 
cal times, except near the outer edge of the disk. Note that 
if hydrostatic equilibrium was not established (i.e. the disk 
was still expanding in the vertical direction, following the 
mass-loss or recoil), the disk would be thinner, and the 3D 
densities would be even higher than our estimates, increas- 
ing the predicted luminosities. Finally, note that under the 
assumption of hydrostatic equilibrium, and with p = £///, 
the Bremsstrahlung integral becomes independent of tem- 
perature, JdVp 2 T 1/2 = const JdrJl 2 r~ 1/2 . 

The solid curve in bottom panel in Figure [14] shows 
the evolution of Z/Brem in the mass-loss only, constant den- 
sity disk runs. The light curve exhibits oscillations within 
the first 100 days after the merger, qualitatively similar to 
the dips and troughs found, on shorter time-scales, by both 
O'Neill et al] (|2009[) and |Megevand et al.| (|2009[). In par- 



ticular, as in the those two studies, the entire disk remains 
dimmer than prior to the mass-loss: the luminosity of the 
disk decays by « 20% over the course of a year. To under- 
stand this result better, we computed the time-dependent 
mean surface density (£} of the simulated disk between 
10 3 — 10 4 7? s . As noted above, the disk as a whole expands af- 
ter the mass-loss, due to the sudden decrease in gravitational 
force from the central point mass. We show, by the dashed 
curve in the bottom panel of Figure [l4| the Bremsstrahlung 
luminosity that would result from a disk with the time- 
dependent surface density (£)t. The departure of the dashed 
curve from unity is therefore due entirely to changes in the 
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Figure 14. Lower panel: Thermal Bremsstrahlung emission from 
an adiabatic, constant-density disk that experiences the effect of 
5% mass loss and no kick (Mad3). The luminosity is normal- 
ized to its steady pre— merger value. Initially, the light curve ex- 
hibits variability, qualitatively similar to the dips and troughs 
found, on shorter time-scales, by both [CVNeill et a.l.| ( |200"9] | and 
|Megevand et al.| p009\ . The dashed curve shows, for compar- 
ison, the Bremsstrahlung luminosity that would result from a 
disk with uniform surface density and the same total mass. The 
prominent dip at t > 100 days is caused by f5i20% of the orig- 
inal disk material exiting the simulation volume. Upper panel: 
Bremsstrahlung emission from our fiducial isothermal disk, expe- 
riencing both mass-loss and a kick (MKiso3). The dashed curve 
shows the corresponding luminosity from a uniform disk, as in the 
lower panel. The disk brightens, mostly due to the sharp density 
contrasts that develop. 



total disk mass within the region 10 3 - W 4 R S , whereas the 
difference between the solid and dashed curve is caused en- 
tirely by inhomogeneities in the gas surface density in this 
region. As these two curves show, at early times, the drop 
in the luminosity is due mostly to inhomogeneities, whereas 
at late times, the dimming is caused by approximately 20% 
of the disk mass exiting the simulated region. Note that at 
late times, the inhomogeneities actually cause an increase 
in the luminosity. At early times, inhomogeneities cause the 
luminosity to decrease. While naively one expects that in- 
homogeneities always boost the luminosity, we note that in 
the adiabatic case, the disk is effectively allowed to expand 
in the vertical direction, and inhomogeneities can cause a 
small overall decrease in the luminosity. 

In comparison, the top panel in Figure \l4\ shows LBrem 
in our fiducial isothermal run (Kiso3), which includes both 
a kick and a mass-loss. As in the bottom panel, the dashed 
curve shows the luminosity from a uniform disk with the 
same total mass. As this figure shows, the luminosity signif- 
icantly increases with time, and this increase is mostly due 
to the inhomogeneities (note that the density contrasts are 
much higher than in the adiabatic case, and the isothermal 
disk does not expand in the vertical direction) . Interestingly, 
in this run that includes the kick, the total mass within the 
10 3 - 10 4 R S region increases, rather than decreases. 

In summary, our constant-density, adiabatic disks, 
which include only mass-loss and no kick, show the same 
qualitative behavior as in |Q'Neillet af]p009] ) and |Megevand| 
et al. (20091, namely oscillations in the luminosity, and 
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a modest overall dimming of the disk, following the BH 
merger. However, using an isothermal, rather than an adia- 
batic equation of state and considering the impact of a kick, 
in addition to the mass-loss, tend to counter this dimming, 
and when combined, produce, instead a significant brighten- 
ing of the disk - due mostly to the stronger density contrasts 
that develop in the isothermal case. 
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4 CONCLUSIONS AND FUTURE 
DIRECTIONS 

In this paper, we used two-dimensional hydrodynamical 
simulations of a circumbinary disk, to follow the effects of a 
velocity recoil and mass-loss of the central black hole binary, 
following the merger of the two black holes. From a suite of 
runs, we are able to draw two basic conclusions. 

First, the outward-propagating spiral shocks that de- 
velop in our simulations follow a pattern very similar to the 



caustics identified in collisionless disk (Lippai et al. 20081 



Gas pressure has a modest overall impact on the propaga- 
tion and on the overall morphology of the shocks. On the 
other hand, we find that the gas pressure, and the assumed 
equation of state, has a strong impact on the overdensities 
that develop and on the strengths of the shocks: isother- 
mal, low-temperature disks have larger overdensities and 
stronger shocks. 

Second, we have estimated an upper limit on the lumi- 
nosity emerging from the disk experiencing a BH kick with 
«kick = 530 km s _1 , by measuring the effective dissipation 
that occurs, implicitly, in the simulations when an isother- 
mal condition is imposed on the gas. The resulting luminos- 
ity is of order 10% of the Eddington limit for the 1O 6 M BH 
used in the simulation, which suggests that the after-glow 
may be bright enough to be detectable. If the pre-merger 
disk has a luminosity below that of a standard steady-state 
thin accretion disk (due to the evacuation of the inner re- 
gions of the disk), then the merger-induced kick will cause 
a significant (order-of-magnitude, or larger) brightening of 
the disk. We also estimated the effective black-body tem- 
perature of the radiation emerging from the optically thick 
disk. We found that as the disk brightens, the characteris- 
tic frequency increases with time, possibly offering a unique 
signature of the kick-induced emission. When the accompa- 
nying mass-loss of the merger remnant is included in our 
simulations, the density contrasts and luminosity from the 
spiral shocks decrease somewhat, but do not change dramat- 
ically in overall behavior. 

While our results are encouraging, and suggest that the 
EM signature of the disturbed circumbinary disk may be 
detectable, this conclusion has to be verified in the future 
in improved models. In particular, a better estimate of the 
thermodynamics of the disk, with a realistic treatment of ra- 
diative cooling, should reveal how close the luminosity is to 
the values we obtained here, using the isothermal assump- 
tion; by using three-dimensional simulations that resolve the 
vertical disk structure, and by following the vertical trans- 
fer of the radiation produced should clarify the robustness 
of our conclusions about the evolution of the characteristic 
emergent frequency. 
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